/**
 * Description:
 * Nearest neighbour functionality using multiple input points
 * 使用多个输入点的最近邻功能
 * (没有测试)
 */

#include "eccodes.h"

static void usage(const char* prog)
{
    fprintf(stderr, "Usage: %s latlon_file grib_orography grib_file grib_file ...\n", prog);
    fprintf(stderr, "       The latlon_file should have 3 columns: num lat lon\n");
    fprintf(stderr, "       The grib_orography file is treated as the land-sea mask\n");
    fprintf(stderr, "\n");
    exit(1);
}

int sfx_grib_nearest_multiple()
{
    int err = 0;

    // Input lat/lon file should have 3 columns:
    //  number   latitude  longitude
    char* fname = "argv[1]";

    int npoints = 0;
    long iid    = 0;
    float lat, lon;

    FILE* fin = fopen(fname, "rt");  // Open in text mode
    while (fscanf(fin, "%ld %g %g", &iid, &lat, &lon) != EOF) {
        npoints++;
    }
    fclose(fin);

    if (npoints == 0) {
        fprintf(stderr, "No input points found!\n");
        exit(1);
    }

    long* id           = (long*)calloc(npoints, sizeof(long));
    double* vlat       = (double*)calloc(npoints, sizeof(double));
    double* vlon       = (double*)calloc(npoints, sizeof(double));
    double* outlats    = (double*)calloc(npoints, sizeof(double));
    double* outlons    = (double*)calloc(npoints, sizeof(double));
    double* values     = (double*)calloc(npoints, sizeof(double));
    double* lsm_values = (double*)calloc(npoints, sizeof(double));
    double* distances  = (double*)calloc(npoints, sizeof(double));
    int* indexes       = (int*)calloc(npoints, sizeof(int));

    //
    fin = fopen(fname, "rt");  // Open in text mode

    int i = 0;
    while (i < npoints && fscanf(fin, "%ld %g %g", &iid, &lat, &lon) != EOF) {
        id[i]   = iid;
        vlat[i] = lat;
        while (lon < 0) {
            lon += 360;
        }
        vlon[i] = lon;
        i++;
    }
    fclose(fin);

    // The first GRIB file on the arguments is treated as the land-sea mask file
    fname = "argv[2]";
    fin   = fopen(fname, "rb");  // Open GRIB in binary mode

    codes_handle* h = codes_handle_new_from_file(0, fin, PRODUCT_GRIB, &err);
    CODES_CHECK(err, 0);

    const int is_lsm = 1;  // ?
    codes_grib_nearest_find_multiple(h, is_lsm, vlat, vlon, npoints,
                                     outlats, outlons, lsm_values, distances, indexes);

    codes_handle_delete(h);
    fclose(fin);

    //
    int argc       = 3;
    size_t len     = 0;
    char* date     = NULL;
    char* time     = NULL;
    long step      = 0;
    long parameter = 0;
    for (int n = 3; n <= argc - 1; n++) {
        //
        fname = "argv[n]";
        fin   = fopen(fname, "rt");
        while ((h = codes_handle_new_from_file(0, fin, PRODUCT_GRIB, &err)) != NULL) {
            codes_get_double_elements(h, "values", indexes, npoints, values);
            //
            CODES_CHECK(codes_get_length(h, "date", &len), 0);
            codes_get_string(h, "date", date, &len);
            CODES_CHECK(codes_get_length(h, "time", &len), 0);
            codes_get_string(h, "time", time, &len);
            codes_get_long(h, "step", &step);
            codes_get_long(h, "parameter", &parameter);
            printf("# %s %s %ld %ld\n", date, time, step, parameter);

            for (i = 0; i < npoints; i++) {
                printf("%ld %.2f %.2f %.2f %.2f %d\n",
                       id[i], outlats[i], outlons[i],
                       lsm_values[i], values[i], indexes[i]);
            }

            codes_handle_delete(h);
        }

        fclose(fin);
    }

    free(id);
    free(vlat);
    free(vlon);
    free(outlats);
    free(outlons);
    free(values);
    free(lsm_values);
    free(distances);
    free(indexes);

    return 0;
}
